In [1]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mplt
import seaborn as sns
import arviz as az
import datetime
In [2]:
import fit_arima
In [3]:
import open_data
df = open_data.open()
In [4]:
ritorno = fit_arima.compute(2, 1, catene=4, samples_per_chain=2500, burnin=1000)
21:10:52 - cmdstanpy - INFO - CmdStan start processing
chain 1 |          | 00:00 Status
chain 2 |          | 00:00 Status
chain 3 |          | 00:00 Status
chain 4 |          | 00:00 Status
                                                                                                                                                                                                                                                                                                                                
21:38:05 - cmdstanpy - INFO - CmdStan done processing.
21:38:05 - cmdstanpy - WARNING - Non-fatal error during sampling:
Exception: normal_lpdf: Scale parameter is -0.111703, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61)
	Exception: code_model_namespace::log_prob: phi[2][11] is 1, but must be less than or equal to 1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 83, column 2 to column 47)
	Exception: code_model_namespace::log_prob: phi[1][1] is -nan, but must be greater than or equal to -1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 83, column 2 to column 47)
	Exception: code_model_namespace::log_prob: phi[1][33] is 1, but must be less than or equal to 1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 83, column 2 to column 47)
Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 157, column 4 to column 34)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 157, column 4 to column 34)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 157, column 4 to column 34)
	Exception: code_model_namespace::log_prob: phi[2][8] is -nan, but must be greater than or equal to -1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 83, column 2 to column 47)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 157, column 4 to column 34)
	Exception: normal_lpdf: Scale parameter is -0.00468339, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61)
	Exception: normal_lpdf: Scale parameter is -0.00781159, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61)
Exception: normal_lpdf: Scale parameter is -19.3346, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61)
	Exception: normal_lpdf: Scale parameter is -20.4458, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61)
	Exception: normal_lpdf: Scale parameter is -1.84055, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61)
	Exception: code_model_namespace::log_prob: phi[2][4] is -nan, but must be greater than or equal to -1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 83, column 2 to column 47)
	Exception: normal_lpdf: Scale parameter is -9.85741, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61)
	Exception: normal_lpdf: Scale parameter is -0.00793045, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61)
Exception: normal_lpdf: Scale parameter is -1.48033, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61)
	Exception: normal_lpdf: Scale parameter is -0.370182, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61)
	Exception: normal_lpdf: Scale parameter is -3.13815, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61)
	Exception: normal_lpdf: Scale parameter is -2.7382, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61)
	Exception: normal_lpdf: Scale parameter is -0.274406, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61)
	Exception: code_model_namespace::log_prob: phi[1][10] is -nan, but must be greater than or equal to -1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 83, column 2 to column 47)
	Exception: normal_lpdf: Scale parameter is -3.44707, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61)
	Exception: normal_lpdf: Scale parameter is -0.00909615, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61)
	Exception: normal_lpdf: Scale parameter is -0.0018971, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61)
Consider re-running with show_console=True if the above output is unclear!

/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/fit_arima.py:37: FutureWarning: Indexing with a float is deprecated, and will raise an IndexError in pandas 2.0. You can manually convert to an integer key instead.
  y_missing_list.append({'Stazione':df.columns[v[i]],'Data':df.index[u[i]],'Samples':inference_data.posterior.w.values[:,:,i].reshape(catene*samples_per_chain)})

Tempo di computazione: 27:12

In [5]:
res = az.loo(ritorno['inference_data'], var_name="log_lik", pointwise=True, scale='deviance')
res
/home/br1/PythonProjects/PythonStats/.venv/lib/python3.10/site-packages/arviz/stats/stats.py:803: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.7 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(
Out[5]:
Computed from 10000 posterior samples and 17191 observations log-likelihood matrix.

             Estimate       SE
deviance_loo -14506.23   309.87
p_loo          296.12        -

There has been a warning during the calculation. Please check the results.
------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.5]   (good)     17112   99.5%
 (0.5, 0.7]   (ok)          68    0.4%
   (0.7, 1]   (bad)         11    0.1%
   (1, Inf)   (very bad)     0    0.0%
In [6]:
aux_plt = az.plot_khat(res, show_bins=True, figsize=(7,5))
In [7]:
res = az.waic(ritorno['inference_data'], var_name="log_lik", scale='deviance')
res
/home/br1/PythonProjects/PythonStats/.venv/lib/python3.10/site-packages/arviz/stats/stats.py:1645: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. 
See http://arxiv.org/abs/1507.04544 for details
  warnings.warn(
Out[7]:
Computed from 10000 posterior samples and 17191 observations log-likelihood matrix.

              Estimate       SE
deviance_waic -14528.97   309.61
p_waic          284.76        -

There has been a warning during the calculation. Please check the results.
In [8]:
plt_aux = az.plot_trace(ritorno['inference_data'], var_names=['phi','theta','hyper_gamma_phi_m','hyper_gamma_phi_s'], divergences=True, figsize=(15,15))
In [9]:
plt_aux = az.plot_trace(ritorno['inference_data'], var_names=['mh_sigma','sh_sigma','sigma'], divergences=True, figsize=(15,10))
In [10]:
plt_aux = az.plot_trace(ritorno['inference_data'], var_names=['hyper_y_start_m','hyper_y_start_s','c'], divergences=True, figsize=(15,10))
In [11]:
az.summary(ritorno['inference_data'],var_names=['y_start'])
Out[11]:
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
y_start[0, 0] 1.422 0.089 1.251 1.587 0.002 0.001 2935.0 3563.0 1.0
y_start[0, 1] 1.419 0.091 1.248 1.591 0.002 0.001 3125.0 3114.0 1.0
y_start[0, 2] 1.422 0.090 1.252 1.593 0.002 0.001 3105.0 4272.0 1.0
y_start[0, 3] 1.424 0.090 1.253 1.594 0.002 0.001 3232.0 3953.0 1.0
y_start[0, 4] 1.427 0.092 1.245 1.595 0.002 0.001 2847.0 3181.0 1.0
... ... ... ... ... ... ... ... ... ...
y_start[2, 44] 1.445 0.054 1.345 1.548 0.001 0.001 1875.0 2615.0 1.0
y_start[2, 45] 1.442 0.055 1.340 1.549 0.001 0.001 1967.0 2244.0 1.0
y_start[2, 46] 1.438 0.054 1.334 1.535 0.001 0.001 2333.0 3204.0 1.0
y_start[2, 47] 1.452 0.055 1.348 1.556 0.001 0.001 1877.0 2748.0 1.0
y_start[2, 48] 1.459 0.057 1.360 1.572 0.001 0.001 1908.0 2766.0 1.0

147 rows × 9 columns

In [12]:
az.summary(ritorno['inference_data'], var_names=['hyper_y_start_m', 'hyper_y_start_s'])
Out[12]:
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
hyper_y_start_m 1.432 0.038 1.364 1.506 0.002 0.001 578.0 852.0 1.00
hyper_y_start_s 0.082 0.012 0.061 0.106 0.000 0.000 688.0 1692.0 1.01
In [13]:
aux_shape = ritorno['inference_data'].posterior.gamma_th.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.gamma_th.values.reshape((aux_shape[0]*aux_shape[1])), legend=False)
Out[13]:
<AxesSubplot: ylabel='Density'>
In [14]:
aux_shape = ritorno['inference_data'].posterior.phi.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.phi.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2],aux_shape[3]))[:,0,:], legend=False)
Out[14]:
<AxesSubplot: ylabel='Density'>
In [15]:
aux_shape = ritorno['inference_data'].posterior.phi.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.phi.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2],aux_shape[3]))[:,1,:], legend=False)
Out[15]:
<AxesSubplot: ylabel='Density'>
In [16]:
aux_shape = ritorno['inference_data'].posterior.c.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.c.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)
Out[16]:
<AxesSubplot: ylabel='Density'>
In [17]:
aux_shape = ritorno['inference_data'].posterior.annual_mean.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.annual_mean.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)
Out[17]:
<AxesSubplot: ylabel='Density'>
In [18]:
aux_shape = ritorno['inference_data'].posterior.annual_median.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.annual_median.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)
Out[18]:
<AxesSubplot: ylabel='Density'>
In [19]:
aux_shape = ritorno['inference_data'].posterior.annual_max.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.annual_max.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)
Out[19]:
<AxesSubplot: ylabel='Density'>
In [20]:
valori_31_dic = []
for dizionario in ritorno['missing']:
    if dizionario['Data'] == '2018-12-31':
        valori_31_dic.append(dizionario['Samples'])
        
res_plot = sns.kdeplot(valori_31_dic, legend=False)
In [21]:
aux_shape = ritorno['inference_data'].posterior.annual_days_over_threshold.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.annual_days_over_threshold.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)
Out[21]:
<AxesSubplot: ylabel='Density'>
In [22]:
pd.set_option('display.max_columns', 500)
aux_results = pd.DataFrame(ritorno['inference_data'].posterior.annual_days_over_threshold.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), columns=df.columns).describe()
aux_results.sort_values('mean', axis=1, ascending=True)
Out[22]:
CASTELLUCCIO CORTE BRUGNATELLA SAVIGNANO DI RIGO FEBBIO SAN LEO VERUCCHIO BADIA GIARDINI MARGHERITA MARECCHIA SAN PIETRO CAPOFIUME PARCO BERTOZZI DELTA CERVIA LUGAGNANO DE AMICIS SAN LAZZARO PARCO RESISTENZA FRANCHINI-ANGELONI GHERARDI VIA CHIARINI GAVELLO SAVIGNANO BESENZONE ROMA VILLA FULVIA PORTA SAN FELICE CAORLE CASTELLARANO CENTO ZALAMELLA PARCO MONTECUCCO PARCO EDILCARANI BOGOLESE SARAGAT PARADIGNA MALCANTONE S. LAZZARO REMESINA S. ROCCO CENO ISONZO PARCO FERRARI GIORDANI-FARNESE CITTADELLA SAN FRANCESCO FLAMINIA MONTEBELLO GIARDINI GERBIDO TIMAVO
count 10000.00000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.00000 10000.000000 10000.000000 10000.000000 10000.00000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.00000 10000.00000 10000.00000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.00000 10000.000000 10000.000000 10000.000000
mean 0.10200 0.115000 0.118500 0.302400 1.567600 4.877900 7.790400 13.415600 15.608700 15.637100 17.239500 17.616200 17.707000 17.769500 18.701500 19.475400 21.104400 21.166500 21.361500 23.260100 23.932500 23.949300 28.18090 28.534400 28.738900 29.321800 29.62880 30.606900 30.616000 30.740400 31.929900 32.091500 34.509200 35.694200 35.87640 36.64660 36.97600 38.350700 41.990600 42.731600 42.800400 46.467400 46.507200 50.969500 51.176500 54.85700 60.543300 63.666400 70.191400
std 0.36443 0.394703 0.398087 0.680147 1.962961 4.054311 5.462239 8.325383 9.067712 9.744521 10.226974 9.588717 9.695842 9.423346 10.193998 10.424829 11.205492 11.605751 11.138134 11.271324 11.896567 12.455571 13.45773 13.180083 13.040335 14.381499 12.01704 14.105356 14.220372 14.289571 12.852171 13.487006 14.189815 15.221756 15.03205 15.35704 14.94129 15.436954 16.984556 16.497745 16.676348 17.410958 17.140721 17.009283 18.785703 18.80313 19.254184 22.584466 21.885408
min 0.00000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.00000 0.000000 0.000000 0.000000 2.00000 1.000000 1.000000 0.000000 2.000000 1.000000 1.000000 2.000000 0.00000 2.00000 2.00000 3.000000 4.000000 1.000000 1.000000 4.000000 7.000000 9.000000 8.000000 8.00000 10.000000 5.000000 9.000000
25% 0.00000 0.000000 0.000000 0.000000 0.000000 2.000000 4.000000 7.000000 9.000000 9.000000 10.000000 11.000000 11.000000 11.000000 11.000000 12.000000 13.000000 13.000000 13.000000 15.000000 15.000000 15.000000 19.00000 19.000000 19.000000 19.000000 21.00000 20.000000 20.000000 20.000000 23.000000 22.000000 24.000000 25.000000 25.00000 26.00000 26.00000 27.000000 30.000000 31.000000 31.000000 34.000000 34.000000 39.000000 38.000000 41.00000 47.000000 48.000000 55.000000
50% 0.00000 0.000000 0.000000 0.000000 1.000000 4.000000 7.000000 12.000000 14.000000 14.000000 15.000000 16.000000 16.000000 16.000000 17.000000 18.000000 19.000000 19.000000 20.000000 21.000000 22.000000 22.000000 26.00000 27.000000 27.000000 27.000000 28.00000 29.000000 28.000000 29.000000 30.000000 31.000000 33.000000 34.000000 34.00000 35.00000 35.00000 36.000000 40.000000 41.000000 41.000000 45.000000 45.000000 50.000000 49.000000 53.00000 59.000000 61.000000 69.000000
75% 0.00000 0.000000 0.000000 0.000000 2.000000 7.000000 10.000000 18.000000 20.250000 21.000000 22.000000 23.000000 23.000000 23.000000 24.000000 25.000000 27.000000 27.000000 27.000000 30.000000 31.000000 31.000000 35.00000 36.000000 36.000000 37.000000 37.00000 38.250000 39.000000 39.000000 39.000000 40.000000 43.000000 44.000000 44.25000 45.00000 46.00000 47.000000 52.000000 52.000000 53.000000 57.000000 57.000000 61.000000 62.000000 66.00000 72.000000 77.000000 83.000000
max 4.00000 7.000000 6.000000 8.000000 21.000000 38.000000 41.000000 80.000000 84.000000 92.000000 91.000000 83.000000 86.000000 85.000000 83.000000 104.000000 102.000000 112.000000 85.000000 101.000000 137.000000 105.000000 130.00000 166.000000 101.000000 111.000000 98.00000 121.000000 119.000000 124.000000 100.000000 119.000000 120.000000 130.000000 130.00000 144.00000 126.00000 120.000000 153.000000 156.000000 141.000000 155.000000 157.000000 139.000000 159.000000 156.00000 163.000000 199.000000 218.000000
In [23]:
pd.set_option('display.max_columns', 500)
aux_results = pd.DataFrame(ritorno['inference_data'].posterior.is_over_daily_limit.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), columns=df.columns).describe()
aux_results.sort_values('mean', axis=1, ascending=True)
Out[23]:
CASTELLUCCIO FEBBIO CORTE BRUGNATELLA SAVIGNANO DI RIGO SAN LEO VERUCCHIO BADIA GIARDINI MARGHERITA MARECCHIA SAN PIETRO CAPOFIUME DE AMICIS LUGAGNANO DELTA CERVIA PARCO BERTOZZI SAN LAZZARO PARCO RESISTENZA FRANCHINI-ANGELONI GHERARDI VIA CHIARINI GAVELLO SAVIGNANO BESENZONE ROMA VILLA FULVIA PORTA SAN FELICE CASTELLARANO CAORLE CENTO PARCO MONTECUCCO ZALAMELLA PARCO EDILCARANI BOGOLESE SARAGAT PARADIGNA MALCANTONE S. LAZZARO REMESINA S. ROCCO CENO ISONZO PARCO FERRARI GIORDANI-FARNESE CITTADELLA FLAMINIA SAN FRANCESCO MONTEBELLO GERBIDO GIARDINI TIMAVO
count 10000.0 10000.0 10000.0 10000.0 10000.0 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.00000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.00000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.00000 10000.000000 10000.000000 10000.000000 10000.000000 10000.00000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000
mean 0.0 0.0 0.0 0.0 0.0 0.000200 0.001300 0.019700 0.034300 0.041900 0.045500 0.051700 0.05210 0.057100 0.066200 0.078100 0.105000 0.108900 0.111100 0.135500 0.157300 0.157400 0.249100 0.258900 0.26580 0.275400 0.287700 0.313500 0.316800 0.319500 0.34190 0.356600 0.426600 0.452100 0.453500 0.47890 0.485400 0.524300 0.610400 0.635600 0.638300 0.717900 0.721100 0.794600 0.816300 0.858800 0.912300 0.916200 0.961700
std 0.0 0.0 0.0 0.0 0.0 0.014141 0.036034 0.138974 0.182008 0.200371 0.208408 0.221432 0.22224 0.232045 0.248644 0.268342 0.306569 0.311529 0.314272 0.342274 0.364102 0.364196 0.432513 0.438052 0.44178 0.446738 0.452713 0.463939 0.465252 0.466306 0.47437 0.479019 0.494608 0.497725 0.497858 0.49958 0.499812 0.499434 0.487684 0.481286 0.480517 0.450044 0.448481 0.404014 0.387259 0.348245 0.282872 0.277101 0.191929
min 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
25% 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000
50% 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000
75% 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000 1.00000 1.000000 1.000000 1.000000 1.000000 1.000000 1.00000 1.000000 1.000000 1.000000 1.000000 1.00000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000
max 0.0 0.0 0.0 0.0 0.0 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.00000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.00000 1.000000 1.000000 1.000000 1.000000 1.000000 1.00000 1.000000 1.000000 1.000000 1.000000 1.00000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000
In [24]:
aux_results.loc['mean',:].plot(kind='barh', figsize=(15,15), xticks=np.linspace(0,1,11))
#res = plt.xticks(rotation = 30)
Out[24]:
<AxesSubplot: >
In [25]:
pd.set_option('display.max_columns', 500)
aux_results = pd.DataFrame(ritorno['inference_data'].posterior.is_over_annual_limit.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), columns=df.columns).describe()
aux_results.sort_values('mean', axis=1, ascending=True)
Out[25]:
CASTELLUCCIO ROMA CASTELLARANO PARCO EDILCARANI DELTA CERVIA BESENZONE PORTA SAN FELICE GAVELLO ZALAMELLA PARCO MONTECUCCO CENTO BOGOLESE SARAGAT REMESINA S. ROCCO GHERARDI LUGAGNANO CAORLE SAN PIETRO CAPOFIUME FEBBIO CORTE BRUGNATELLA SAN LAZZARO SAN LEO VERUCCHIO BADIA GIARDINI MARGHERITA SAVIGNANO DI RIGO PARCO BERTOZZI MARECCHIA DE AMICIS VIA CHIARINI FRANCHINI-ANGELONI PARCO RESISTENZA VILLA FULVIA MALCANTONE SAVIGNANO S. LAZZARO CENO ISONZO PARADIGNA SAN FRANCESCO PARCO FERRARI GIORDANI-FARNESE CITTADELLA FLAMINIA MONTEBELLO GIARDINI GERBIDO TIMAVO
count 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0 10000.0000 10000.0000 10000.0000 10000.0000 10000.0000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000
mean 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0001 0.0001 0.0001 0.0001 0.0001 0.000200 0.000200 0.000300 0.000400 0.000400 0.000600 0.001400 0.001700 0.002400 0.014100 0.014600
std 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0100 0.0100 0.0100 0.0100 0.0100 0.014141 0.014141 0.017319 0.019997 0.019997 0.024489 0.037392 0.041198 0.048933 0.117909 0.119951
min 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0000 0.0000 0.0000 0.0000 0.0000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
25% 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0000 0.0000 0.0000 0.0000 0.0000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
50% 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0000 0.0000 0.0000 0.0000 0.0000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
75% 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0000 0.0000 0.0000 0.0000 0.0000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
max 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0000 1.0000 1.0000 1.0000 1.0000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000
In [26]:
aux_results.loc['mean',:].plot(kind='barh', figsize=(15,15), xticks=np.linspace(0,1,11))
#res = plt.xticks(rotation = 30)
Out[26]:
<AxesSubplot: >
In [27]:
y_post_pred = ritorno['inference_data'].posterior.y_post_pred.to_numpy()
y_post_pred = y_post_pred.reshape(y_post_pred.shape[0]*y_post_pred.shape[1], y_post_pred.shape[2], y_post_pred.shape[3])
y_post_pred.shape
Out[27]:
(10000, 365, 49)
In [28]:
import ipywidgets as widgets
from ipywidgets import HBox, VBox
from IPython.display import display
In [29]:
@widgets.interact(stazione = df.columns)
def f(stazione):
    col_map = sns.light_palette((20, 75, 70), input='husl', as_cmap=True)

    plt.figure(figsize=(14,8))
    ax = plt.subplot(1,1,1)

    idx_stazione = df.columns.to_list().index(stazione)
    
    """
    sns.lineplot(np.transpose(y_post_pred[0:50,:,idx_stazione]))
    """

    lower_lim = np.zeros(len(df.index))
    upper_lim = np.zeros(len(df.index))    
    for i in range(len(df.index)):
        lower_lim[i] = np.percentile(y_post_pred[:,i,idx_stazione],2.5)
        upper_lim[i] = np.percentile(y_post_pred[:,i,idx_stazione],97.5)
    
    
#    sns.lineplot(lower_lim)
#    sns.lineplot(upper_lim)
    
    for i in range(len(df.index)):
        ax.add_patch(mplt.patches.Rectangle((i,lower_lim[i]),1,upper_lim[i]-lower_lim[i], fill=True, color=col_map(180)))

    sns.lineplot(df[stazione])    
    sns.lineplot(np.mean(y_post_pred[:,:,idx_stazione], axis=0))


    index = 0
    for line in ax.get_lines():
        if index == 0:
            col_map = sns.dark_palette((230,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(155))
        else:
            col_map = sns.dark_palette((10,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(175))
        index += 1
        
    plt.ylim(bottom=0,top=2.5)

    primo_giorno = datetime.date(2018,1,1)
    date_da_segnare = []
    date_da_segnare_posizioni = []

    for i in range(12):
        date_da_segnare.append(datetime.date(2018,i+1,1))
        date_da_segnare_posizioni.append((date_da_segnare[2*i] - primo_giorno).days)
        date_da_segnare.append(datetime.date(2018,i+1,15))
        date_da_segnare_posizioni.append((date_da_segnare[2*i+1] - primo_giorno).days)
        date_da_segnare[2*i] = date_da_segnare[2*i].isoformat()
        date_da_segnare[2*i+1] = date_da_segnare[2*i+1].isoformat()

    plt.xticks(date_da_segnare_posizioni,date_da_segnare,rotation=65)
    plt.tick_params(
        axis='x',          # changes apply to the x-axis
        which='both',      # both major and minor ticks are affected
        bottom=True,      # ticks along the bottom edge are off
        top=False,         # ticks along the top edge are off
        labelbottom=True)
    plt.grid()
    
    """
    ax.get_legend().remove()
    col_vals = np.linspace(1,255,num=len(df.columns))
    index = 0
    for line in ax.get_lines():
        if(index == len(ax.get_lines()) - 1):
            col_map = sns.dark_palette((230,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(175))
            continue
        if(index == len(ax.get_lines()) - 2):
            col_map = sns.dark_palette((120,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(175))
            index += 1
            continue

        #line.set_c(col_map(int(np.round(col_vals[index]))))
        line.set_c(col_map(220))
        line.set_alpha(0.3)
        index += 1
    """
    
    

    plt.show()
interactive(children=(Dropdown(description='stazione', options=('CASTELLUCCIO', 'FEBBIO', 'CORTE BRUGNATELLA',…
In [30]:
def f(stazione):
    col_map = sns.light_palette((20, 75, 70), input='husl', as_cmap=True)

    plt.figure(figsize=(14,8))
    ax = plt.subplot(1,1,1)

    idx_stazione = df.columns.to_list().index(stazione)
    
    """
    sns.lineplot(np.transpose(y_post_pred[0:50,:,idx_stazione]))
    """

    lower_lim = np.zeros(len(df.index))
    upper_lim = np.zeros(len(df.index))    
    for i in range(len(df.index)):
        lower_lim[i] = np.percentile(y_post_pred[:,i,idx_stazione],2.5)
        upper_lim[i] = np.percentile(y_post_pred[:,i,idx_stazione],97.5)
    
    
#    sns.lineplot(lower_lim)
#    sns.lineplot(upper_lim)
    
    for i in range(len(df.index)):
        ax.add_patch(mplt.patches.Rectangle((i,lower_lim[i]),1,upper_lim[i]-lower_lim[i], fill=True, color=col_map(180)))

    sns.lineplot(df[stazione])    
    sns.lineplot(np.mean(y_post_pred[:,:,idx_stazione], axis=0))


    index = 0
    for line in ax.get_lines():
        if index == 0:
            col_map = sns.dark_palette((230,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(155))
        else:
            col_map = sns.dark_palette((10,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(175))
        index += 1
        
    plt.ylim(bottom=0,top=2.5)

    primo_giorno = datetime.date(2018,1,1)
    date_da_segnare = []
    date_da_segnare_posizioni = []

    for i in range(12):
        date_da_segnare.append(datetime.date(2018,i+1,1))
        date_da_segnare_posizioni.append((date_da_segnare[2*i] - primo_giorno).days)
        date_da_segnare.append(datetime.date(2018,i+1,15))
        date_da_segnare_posizioni.append((date_da_segnare[2*i+1] - primo_giorno).days)
        date_da_segnare[2*i] = date_da_segnare[2*i].isoformat()
        date_da_segnare[2*i+1] = date_da_segnare[2*i+1].isoformat()

    plt.xticks(date_da_segnare_posizioni,date_da_segnare,rotation=65)
    plt.tick_params(
        axis='x',          # changes apply to the x-axis
        which='both',      # both major and minor ticks are affected
        bottom=True,      # ticks along the bottom edge are off
        top=False,         # ticks along the top edge are off
        labelbottom=True)
    plt.grid()
    
    """
    ax.get_legend().remove()
    col_vals = np.linspace(1,255,num=len(df.columns))
    index = 0
    for line in ax.get_lines():
        if(index == len(ax.get_lines()) - 1):
            col_map = sns.dark_palette((230,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(175))
            continue
        if(index == len(ax.get_lines()) - 2):
            col_map = sns.dark_palette((120,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(175))
            index += 1
            continue

        #line.set_c(col_map(int(np.round(col_vals[index]))))
        line.set_c(col_map(220))
        line.set_alpha(0.3)
        index += 1
    """
    
    

    plt.show()
In [31]:
for stazione in df.columns:
    f(stazione)
In [ ]:
 
In [ ]:
 
In [ ]: